Search for: All records

Abstract. Many geodynamical models are formulated in terms of the Stokes equations that are then coupled to other equations. For the numerical solution of the Stokes equations, geodynamics codes over the past decades have used essentially every finite element that has ever been proposed for the solution of this equation, on both triangular/tetrahedral (“simplex”) and quadrilaterals/hexahedral (“hypercube”) meshes. However, in many and perhaps most cases, the specific choice of element does not seem to have been the result of careful benchmarking efforts but based on implementation efficiency or the implementers' background. In a first part of this paper (Thieulot and Bangerth, 2022), we have provided a comprehensive comparison of the accuracy and efficiency of the most widely used hypercube elements for the Stokes equations. We have done so using a number of benchmarks that illustrate “typical” geodynamic situations, specifically taking into account spatially variable viscosities. Our findings there showed that only Taylor–Hood-type elements with either continuous (Q2×Q1) or discontinuous (Q2×P-1) pressure are able to adequately and efficiently approximate the solution of the Stokes equations. In this, the second part of this work, we extend the comparison to simplex meshes. In particular, we compare triangular Taylor–Hood elements against the MINI element and one often referred to as the “Crouzeix–Raviart” element. We compare these choices against the accuracy obtained on hypercube Taylor–Hood elements with approximately the same computational cost. Our results show that, like on hypercubes, the Taylor–Hood element is substantially more accurate and efficient than the other choices. Our results also indicate that hypercube meshes yield slightly more accurate results than simplex meshes but that the difference is relatively small and likely unimportant given that hypercube meshes often lead to slightly denser (and consequently more expensive) matrices. 

This paper provides an overview of the new features of the finite element library deal.II, version 9.6. 

Bayesian methods have been widely used in the last two decades to infer statistical proper- ties of spatially variable coefficients in partial differential equations from measurements of the solutions of these equations. Yet, in many cases the number of variables used to param- eterize these coefficients is large, and oobtaining meaningful statistics of their probability distributions is difficult using simple sampling methods such as the basic Metropolis– Hastings algorithm—in particular, if the inverse problem is ill-conditioned or ill-posed. As a consequence, many advanced sampling methods have been described in the literature that converge faster than Metropolis–Hastings, for example, by exploiting hierarchies of statistical models or hierarchies of discretizations of the underlying differential equation. At the same time, it remains difficult for the reader of the literature to quantify the advantages of these algorithms because there is no commonly used benchmark. This paper presents a benchmark Bayesian inverse problem—namely, the determination of a spatially variable coefficient, discretized by 64 values, in a Poisson equation, based on point mea- surements of the solution—that fills the gap between widely used simple test cases (such as superpositions of Gaussians) and real applications that are difficult to replicate for de- velopers of sampling algorithms. We provide a complete description of the test case and provide an open-source implementation that can serve as the basis for further experiments. We have also computed 2 × 10^11 samples, at a cost of some 30 CPU years, of the poste- rior probability distribution from which we have generated detailed and accurate statistics against which other sampling algorithms can be tested. 

Thehp-adaptive finite element method—where one independently chooses the mesh size (h) and polynomial degree (p) to be used on each cell—has long been known to have better theoretical convergence properties than eitherh- orp-adaptive methods alone. However, it is not widely used, owing at least in part to the difficulty of the underlying algorithms and the lack of widely usable implementations. This is particularly true when used with continuous finite elements. Herein, we discuss algorithms that are necessary for a comprehensive and generic implementation ofhp-adaptive finite element methods on distributed-memory, parallel machines. In particular, we will present a multistage algorithm for the unique enumeration of degrees of freedom suitable for continuous finite element spaces, describe considerations for weighted load balancing, and discuss the transfer of variable size data between processes. We illustrate the performance of our algorithms with numerical examples and demonstrate that they scale reasonably up to at least 16,384 message passage interface processes. We provide a reference implementation of our algorithms as part of the open source librarydeal.II. 

Abstract. Numerical models are a powerful tool for investigating the dynamic processes in the interior of the Earth and other planets, but the reliability and predictive power of these discretized models depends on the numerical method as well as an accurate representation of material properties in space and time. In the specific context of geodynamic models, particle methods have been applied extensively because of their suitability for advection-dominated processes and have been used in applications such as tracking the composition of solid rock and melt in the Earth's mantle, fluids in lithospheric- and crustal-scale models, light elements in the liquid core, and deformation properties like accumulated finite strain or mineral grain size, along with many applications outside the Earth sciences. There have been significant benchmarking efforts to measure the accuracy and convergence behavior of particle methods, but these efforts have largely been limited to instantaneous solutions, or time-dependent models without analytical solutions. As a consequence, there is little understanding about the interplay of particle advection errors and errors introduced in the solution of the underlying transient, nonlinear flow equations. To address these limitations, we present two new dynamic benchmarks for transient Stokes flow with analytical solutions that allow us to quantify the accuracy of various advection methods in nonlinear flow. We use these benchmarks to measure the accuracy of our particle algorithm as implemented in the ASPECT geodynamic modeling software against commonly employed field methods and analytical solutions. In particular, we quantify if an algorithm that is higher-order accurate in time will allow for better overall model accuracy and verify that our algorithm reaches its intended optimal convergence rate. We then document that the observed increased accuracy of higher-order algorithms matters for geodynamic applications with an example of modeling small-scale convection underneath an oceanic plate and show that the predicted place and time of onset of small-scale convection depends significantly on the chosen particle advection method. Descriptions and implementations of our benchmarks are openly available and can be used to verify other advection algorithms. The availability of accurate, scalable, and efficient particle methods as part of the widely used open-source code ASPECT will allow geodynamicists to investigate complex time-dependent geodynamic processes such as elastic deformation, anisotropic fabric development, melt generation and migration, and grain damage.